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PREFACE 


This  work  is  related  to  a  research  program  on  ocean  bottom  characterization  being  conducted  by 
the  Acoustics  Media  Characterization  Branch  at  NRL.  The  program  is  motivated  by  the  recent 
development  of  a  variety  of  bathymetric  measurement  systems  that  resolve  aspects  of  bottom  variabil¬ 
ity  heretofore  beyond  the  realm  of  quantitative  study.  The  general  aim  of  the  research  is  to  develop 
statistical/stochastic  representations  of  ocean  bottom  variability  at  scales  that  are  relevant  to  underwa¬ 
ter  acoustic  systems.  Such  representations  will  allow  the  prediction  of  acoustic  system  performance 
and  lead  to  the  better  design  and  use  of  underwater  acoustic  systems  that  either  utilize  or  are  sensitive 
to  acoustic  interaction  with  the  ocean  bottom. 


IV 


A  METHOD  FOR  VALIDATING  MULTIDIMENSIONAL 
FAST  FOURIER  TRANSFORM  (FFT)  ALGORITHMS 


INTRODUCTION 

In  the  physical  sciences,  many  problems  exist  that  lead  to  the  analysis  of  a  physical  process  in 
terms  of  the  concept  of  a  random  process  (Yaglom  1962).  A  key  aspect  of  the  theory  of  random 
processes  is  the  spectrum  of  the  process,  which  describes  the  distribution  of  energy  associated  with 
the  process  over  wave  number  or  frequency.  The  computation  of  the  spectrum  is  often  carried  out  by 
evaluating  a  discrete  Fourier  transform  using  a  fast  Fourier  transform  (FFT)  algorithm.  This  report 
deals  with  some  problems  that  arise  in  connection  with  FFT  algorithms. 

Many  potential  difficulties  arise  in  connection  with  the  use  of  an  FFT.  1.  The  usual  situation  is 
that  the  algorithm  already  exists  as  an  available  subroutine  on  the  computer  facility;  however,  the 
documentation  of  the  algorithm  may  be  incomplete  and  there  may  be  questions  as  to  which  of  several 
possible  definitions  of  the  transform  is  being  utilized.  2.  Another  situation  that  arises  is  the  writing  of 
an  FFT  algorithm  for  a  computer  facility,  possibly  to  take  advantage  of  a  new  technique.  In  this 
case,  the  specific  algorithm  for  the  transform  is  presumably  known,  but  it  is  desirable  to  have  a 
means  to  check  the  validity  of  the  calculation  once  the  code  has  been  implemented.  3.  A  further 
situation  that  may  arise  is  that  an  FFT  may  be  added  to  a  computer  facility  as  a  built-in  feature  of  an 
array  processor.  Again  the  question  arises  as  to  whether  the  array  processor  is  doing  the  calculation 
correctly. 

We  have  mentioned  three  situations  in  which  it  is  desirable  to  have  a  means  to  check  the  validity 
of  the  calculation  of  the  FFT.  This  is  especially  important  because  such  an  algorthm  will  receive 
widespread  application.  This  report  provides  a  method  for  validating  FFT  algorithms.  The  method  is 
specifically  dealt  with  for  one-dimensional  and  two-dimensional  discrete  Fourier  transforms.  The 
two-dimensional  case  demonstrates  how  the  method  can  be  easily  generalized  to  multidimensional 
transforms  of  any  dimension. 

TERMINOLOGY  AND  BASIC  RELATIONS 

The  discrete  Fourier  transform  is  at  the  core  of  any  procedure  that  evaluates  a  spectrum  from 
digital  data.  In  this  section  we  summarize  those  relationships  that  will  be  needed  in  the  following 
analysis.  A  more  extensive  derivation  and  discussion  of  discrete  Fourier  transform  properties  can -be 
found  in  several  texts  (Brigham  1974,  Oppenheim  and  Schafer  1975). 

For  one-dimensional  data  X(k),  the  discrete  Fourier  transform  is  defined  as 

M-l 

A{m)  =  £  X(l)-2nM/M  (m  =  0, 1, . . .  ,M  —  1) ,  (1) 

k  =0 

where  it  is  assumed  that  the  original  data  consists  of  M  discrete  values.  It  is  useful  to  extend  the 
definition  of  A(m)  and  X (k )  to  all  discrete  values  of  m  and  £  by  a  periodic  extension  of  the  previous 
set  of  values.  The  discrete  transform  A(m)  is  in  general  complex.  For  real  data,  however,  we  have 

X*(k)  =  X(k),  (2) 
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and  this  leads  to 

A{m)  =A*(—m )  =  A*(M  —m), 
for  m  =  0, 1, . . .  ,M  -  1 . 


(3) 


A  basic  relation  that  exists  between  the  transform  and  the  data  is  given  by  Parseval’s  relation, 
which  states  that 

M  — 1  1  M- 1 

£  \x(k)\2  =  T7  £  \A(m)\2.  (4) 

k  =0  M  m  =0 

This  relationship  can  be  thought  of  as  an  energy  relation  and  is  useful  for  establishing  correspon¬ 
dences  between  the  spectrum  of  the  continuous  process  and  the  discrete  Fourier  transform,  which  is 
calculated  from  equally  spaced  measurements  of  the  process. 


We  now  consider  the  two-dimensional  case  where  we  have  data  X(k,l )  given  over  a  two- 
dimensional  grid.  The  discrete  Fourier  transform  is  now  given  by 

M- 1  N-\  „  .  ,  „  „  . 

A(m,n)  =  £  J2x(k’l>>e  e  (5) 

k=  0  1=0 


where  the  data  X(k,l )  is  assumed  given  over  a  set  of  data  points  identified  by  k  =  0, 1,. . .  ,M-1 
and  /  =  0, 1, . . .  ,N  -  1.  It  is  again  useful  to  regard  the  transform  as  defined  for  all  integer  values  of 
m,  n.  When  the  data  X(k,l)  are  real,  we  obtain  the  symmetry  properties 

A*(—m,  — n )  =  A(m  ,n)  =  A*(M  —  m  ,N  —  n) .  (6) 

Parseval’s  relation  for  discrete  data  over  a  two-dimensional  region  takes  the  form 

M- 1  N- 1  1  M-l  N- 1 

£  £  W'".")!  =^77  £  £  I2.  (7) 

m=0  n=0  MiV  k=  0  /=0 


An  important  view  to  take  toward  discrete  Fourier  transforms  of  data  that  are  defined  over  a 
multidimensional  array  of  points  is  to  regard  the  transform  as  a  series  of  one-dimensional  transforms. 
For  example,  the  two-dimensional  transform  defined  in  Eq.  (5)  can  be  written  as 

N- 1  M-i  „  ,  w  „  .  I/xr 

A(m,n)  =  £  {  £  X(k ,l)e~2lTimk/M  }  e~2Tml/N  .  (8) 

1=0  k =0 

In  this  form,  it  can  be  regarded  as  the  result  of  several  one-dimensional  transforms  carried  out  first  on 
the  rows  of  the  data  array  X(k,l).  A  series  of  one-dimensional  transforms  is  then  performed  on  each 
of  the  columns  of  the  intermediate  array  to  yield  the  two-dimensional  transform  of  the  original  array 
of  data.  One  could  just  as  well  begin  with  column  transforms  rather  than  row  transforms  in  Eq.  (8). 
This  point  of  view  is  obviously  capable  of  generalization  to  higher  dimensions  and  is  also  the  basis 
for  computational  algorithms  (Robinson  and  Silvia,  1979).  Thus  one  can  regard  the  one-dimensional 
discrete  Fourier  transform  as  the  fundamental  ingredient  for  the  higher  dimensional  transforms. 

We  consider  one-dimensional  and  two-dimensional  transforms  in  this  report  because  there  are 
some  aspects  of  higher  dimensional  transforms  that  become  clear  only  when  one  deals  with  a 
transform  of  data  defined  over  more  than  one  dimension. 

DESCRIPTION  OF  THE  METHOD 

The  essential  aspect  of  the  method  for  validating  discrete  Fourier  transforms  is  to  use  simple 
discrete  functions  for  which  the  transform  can  be  worked  out  analytically— by  this  we  mean  the  sum¬ 
mation  appearing  in  the  transform  relation  can  be  carried  out  to  obtain  an  expression  involving  ele¬ 
mentary  functions  whose  number  does  not  vary  with  the  number  of  input  values.  One  then  has  a 
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means  for  checking  the  output  of  the  transform  algorithm.  This  procedure  allows  the  transform  algo¬ 
rithm  to  be  viewed  as  a  black  box ,  which  may  be  advantageous  in  some  situations.  We  shall  consider 
several  simple  functions  that  can  be  used  as  input  and  discuss  their  advantages  and  disadvantages  as  a 
test  of  the  transform  algorithm.  The  one-dimensional  and  two-dimensional  transforms  are  considered 
because  there  are  some  potential  difficulties  with  the  multidimensional  transform  that  do  not  arise  with 
the  one-dimensional  transform.  Further,  the  consideration  of  the  two-dimensional  transform  will 
demonstrate  how  the  method  can  be  generalized  to  multidimensional  transforms  of  any  dimension. 


A  simple  test  is  provided  by  Parseval’s  relation,  which  relates  a  property  of  the  original  function 
to  a  property  of  the  transform.  This  relation  can  be  used  to  check  the  correctness  of  the  transform 
algorithm  by  using  a  simple  function  as  input.  Perhaps  the  simplest  input  function  is 

X(k)  =  8(k)  (k  =  0,1,..., M-l),  (9) 

where 

6(k)  =  1  if  k  =  0  (mod M)  (10) 

=  0  otherwise. 


The  discrete  Fourier  transform  is 


which  becomes 


M  — 1 

A{m)  =  £  8(k)e-2ltimk/M  , 
k  =0 


A(m)  =  1  for  m  =  0, 1,. . .  ,M  —  1 . 


Parseval’s  relation  states  that  we  must  have 


M- 1  i  M-l 

E  I X(k)  1 2  =  "77  E  \Mm)\2. 

k=  0  M  m-0 


The  summation  on  the  right  is 


i  M  —  \ 

mE1  =  1- 

m  m  =0 


which  is  also  what  the  summation  of  the  square  of  the  input  data  values  equals. 


This  is  easily  generalized  to  higher  dimensions  as  is  seen  for  the  two-dimensional  case  where 
one  would  use  X(k,l)  =  8(k)8(l)  as  the  input.  The  discrete  Fourier  transform  is 

M  — 1  N- 1 

A(m,n)  =  E  E  S(k)8(l)e-2*imk/M  e~2*inl/N  , 

k =0  1=0 


which  becomes 


A  (m ,  n )  =  1  for  m  =  0, 1, . . . ,  M  —  1 
n  =  0, 1, . . .  ,  iV  —  1 . 

Again  Parseval’s  relation  requires  that 

M  — 1  N- 1  i  M-l  N- 1 

E  E  l*(M)|2  =  "A;  E  E  \A(m,n)  |2, 

k  =0  1=0  Mly  m=  0  n  =0 

where  it  is  easily  verified  that  each  side  of  the  equality  has  the  value  1 . 

These  results  provide  one  means  of  testing  the  validity  of  the  discrete  Fourier  transform  algo¬ 
rithm  by  indicating  what  results  should  be  obtained  from  an  input  of  a  simple  form.  The  disadvan¬ 
tage  of  this  test  is  that  the  resultant  transforms  have  all  elements  equal  to  one  another.  It  is  conceiv¬ 
able  that  a  coding  error  could  produce  a  constant  entry  in  all  array  locations  that  would  then  yield  a 
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false  indication  of  a  correct  result  for  the  simple  preceding  input  data.  It  therefore  seems  desirable  to 
use  a  more  complicated  input  function  to  obtain  a  transform  with  sufficient  variability  that  it  would  be 
unlikely  to  provide  a  false  indicator. 

It  is  important  to  realize  that  Parseval’s  relation  is  a  limited  test  based  on  a  transform  property 
that  is  not  restricted  to  Fourier  transforms.  Some  further  aspects  of  Parseval’s  relation  that  prevent  it 
from  serving  as  a  detailed  test  are  that  the  relation  involves  magnitudes,  which  do  not  preserve  phase 
relations  between  various  coefficients,  and  the  summation  is  insensitive  to  a  permutation  of  the  order 
of  the  coefficients.  These  features  of  Parseval’s  relation  imply  that  the  relation  is  only  useful  as  a 
global  test.  The  most  stringent  test  of  an  FFT  algorithm  will  require  examination  of  the  specific 
values  of  transform  coefficients  and  should  involve  the  use  of  a  fairly  general  input  function. 

We  now  consider  the  advantages  of  the  linear  function 

X(k)=k  (k  =  0,1 . M  -  1).  (11) 


The  transform  is 

M- 1 

Aim)  =  £  ke~2imklM  (m  =  0, 1, . . .  ,M  -  1)  (12) 

k=  0 

and  can  be  evaluated  explicitly  in  terms  of  simple  functions  as  shown  in  the  appendix  of  this  report. 
The  result  is  that 


A  (0) 


,  M  .  M  . 
A(m)  =  -  —  +  i  —  cot 


M(M  -  1) 
2 


ion 

M 


(m  —  1,2,. . .  ,M  —  1). 


(13) 


It  is  seen  from  the  expressions  in  Eq.  (13)  that  the  various  elements  of  the  transform  are  different, 
thus  the  simple  test  function  provides  a  useful  means  for  validating  one-dimensional  discrete  Fourier 
transforms. 


The  test  function  appearing  in  Eq.  (11)  can  easily  be  generalied  to  higher  dimensions.  We 
explicitly  consider  the  two-dimensional  case  in  which 

k  =  0,. . .  ,M  -  1 
/  =  0, . . .  ,N  -  1 

The  two-dimensional  discrete  Fourier  transform  is  given  by 

-2vi mk/M  „ 


X(k,l )  =  kl 


M-\  N- 1 

A(m,n)  =  £  £  Me 

k=  o  1=0 

m  =  0,. . .  ,M  —  1 
n  =  0, . . .  ,N  -  1 


which  can  be  written  as 


A  (m ,  n )  = 

r  "> 

M- 1  .  ,  w 

j^g—lnimk/M 

N- 1 

s. k=0  , 

J=° 

where 


m  =  0,. . .  ,M  —  1 
n  =  0, . . .  ,N  —l 


(14) 


(15) 
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The  two-dimensional  transform  in  Eq.  (15)  is  obviously  the  product  of  one-dimensional  transforms, 
which  can  be  written  out  explicitly  using  the  results  obtained  in  the  appendix  of  this  report.  The 
result  is 


A  (0, 0)  =  -M(M  -  1  )N(N  -  1); 


A(0,n) 

A(m,0) 


J_ 

2 


1) 

_N_ 

,  -N 

r 

7Ttt 

2 

+ 1—  cot 

2 

~N 

-J 

yW(AT-l) 


“ 

r 

~ 

M 

,  .M 

7 Tin 

2 

+ 1—  cot 

2 

M 

v-  -J 

0 n 

(m 


1.2..  .  .  ,N  —  l) , 

1.2..  . .  ,M  —  1), 


and 


where 


A{m,n)  — 


M  ,  .M 
•  — —  +  i—  cot 
2  2 


- 

irm 

^ _ 

N  ,  .N 
—  + 1—  cot 
2  2 


r 

- 

irm 

( 

^ _ 

(16) 

(17) 

(18) 


(19) 


m  =  1,2,. . .  ,M  —  1 
n  =  1,2,...  ,N  -  1  ' 

v  J 

It  is  useful  to  think  of  the  results  in  Eqs.  (16)  to  (19)  as  providing  various  entries  in  a  matrix  as 
shown  in  Fig.  1. 


The  two-dimensional  test  function,  just  as  in  the  one-dimensinoal  case,  leads  to  a  transform  with 
variable  entries  and  provides  a  useful  means  for  validating  two-dimensional  FFT  algorithms.  The  test 
function  can  easily  be  generalized  to  higher  dimensions  in  which  case  the  multidimensional  discrete 
Fourier  transform  of  the  test  function  is  given  by  the  product  of  a  number  of  one-dimensional 
transforms,  entirely  analogous  to  the  two-dimensional  case  we  have  just  considered. 


Fig.  1  —  Matrix  form  of  the  two-dimensional  transform  A  (m ,  n )  that  is 
explicitly  written  out  in  Eqs.  (16)  to  (19).  The  indexes  m  and  n 
correspond  to  row  and  column  indicators,  respectively.  The  reason  for 
exhibiting  the  transform  in  matrix  form  is  to  emphasize  the  structural 
aspects  of  the  transform. 
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The  function  X(k,l )  =kl  we  have  been  considering  seems  adequate  for  testing  the  validity  of 
the  general  case.  However,  it  has  some  disadvantages  when  M  =  N  because  in  that  case  the 
transform  A{m,n )  becomes  a  symmetric  matrix.  Thus  for  the  case  M  =  N,  the  test  function 
X(k,l)  =  kl  has  properties  that  are  not  shared  in  general  by  two-dimensional  input  functions.  The 
symmetry  of  the  transform  makes  it  impossible  to  test  whether  or  not  rows  and  columns  have  been 
incorrectly  interchanged  in  the  program.  This  becomes  an  important  consideration  in  testing  the  va¬ 
lidity  of  algorithms  that  are  specifically  written  for  multidimensional  input  with  M  =  N .  To  have  a 
means  for  testing  this,  we  need  to  consider  some  other  input  functions  that  do  not  lead  to  symmetric 
transforms. 


Consider  the  following  functions: 

X(k,l)  =  ak  +bl,  (20) 

X(k,l)  =  ak  +  b  ,  and  (21) 

X(k,l)=a+bl.  (22) 

In  each  of  these,  a  and  b  are  constants.  We  shall  deal  with  the  general  case  with  M  *  N  although 
our  main  interest  is  in  the  case  when  M  =  N.  Again,  the  functions  selected  are  of  a  sufficiently  sim¬ 
ple  form  that  their  transforms  may  be  worked  out  analytically. 


The  transform  of  the  discrete  function  given  in  Eq.  (20)  is 

M -1  N-l 

A(m,n)  =  £  £  (ak  +  bl)e~2*,mk/M e-2nnl/N  , 

k =0  / =0 


which  can  be  written  as 

A  (m ,  n  )  =  a 

+  b 


S' 

M- 1  .  ,  , 

jQg  —iTrimk  / M 

rN- 1  ^ 

y’s  £—27 xinllN 

l*=° 

J=° 

M  —  1 

^W-l 

£—27 xirnklM 

£  le  ~2T,ml/N 

II  1 
o 

A  _ 

J-  0 

(23) 


The  two-dimensional  transform  is  thus  given  by  evaluating  a  series  of  one-dimensional  transforms. 
Note  that  we  have  the  simple  result  that 

NJ2  e~2lTinl/N  =  Nb(n ) ,  (24) 

1=0 


where 

8(n)  =  lifn  =0  (modA) 
=  0  otherwise. 


The  expression  in  Eq.  (24)  is  the  one-dimensional  transform  of  the  function  that  is  equal  to  unity. 
The  result  in  Eq.  (24)  allows  Eq.  (23)  to  be  written  as 


M  - 1 

r 

N-\  ,  .  ,  „ 

A(m,n)  =  aN 

y^  ^  —lirimk  / M 

5(n )  +  bM 

£  le~  2*mllN 

o 

;  II 

1 

J=° 

8(m). 


(25) 


The  general  structure  of  this  transform  is  best  thought  of  in  terms  of  a  matrix  where  m  designates 
row  location  and  n  designates  column  location.  Figure  2  shows  the  general  form  of  the  result  where 
nonzero  entries  only  occur  in  the  first  row  or  first  column.  These  entries  are  given  by 


A  (0, 0)  =  aNM(m2  0  +bMN(N2 


A  (m ,  0)  =  aN 


M  ,  .M 
— —  +  i—  cot 
2  2 


s-  ~\ 

— 

7r  m 

_ J 

( m  =  1,2,. . .  ,M  —  1) 


(26) 

(27) 
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and 


A(0,n)  =  bM 

N  N 

f  ~N 

7171 

1 

+ 

o 

o 

r-+ 

P~N 

_ 

L  J 

1, 2, ...  ,1V  —  1) , 


where  use  has  been  made  of  the  results  obtained  in  the  appendix  of  this  report. 


(28) 


Fig.  2  —  Matrix  form  of  the  transform  A(m,n)  of  the  function 
X(k,l)  =  ak  +bl  showing  its  structural  form.  Explicit  row  and 
column  entries  are  given  by  Eqs.  (26)  to  (28). 


When  M  =  N,  the  only  difference  between  the  first  row  and  first  column  are  in  the  constants  a 
and  b .  Thus,  while  the  given  test  function  allows  a  test  of  row  and  column  placement,  it  does  not 
involve  much  variability  in  the  entries. 


The  two-dimensional  discrete  Fourier  transforms  of  the  functions  in  Eqs.  (21)  and  (22)  can  be 
evaluated  by  using  the  results  that  were  used  in  the  evaluation  of  the  transform  of  Eq.  (20).  The 
transform  of  Eq.  (21)  leads  to 


A(m  ,n)  =  aN 


rM- 1 

E  ke 


— 2tt  imk  /M 


V- 


k  =0 


8(n )  +  bMN 8(m)8(n) , 


(29) 


while  the  transform  of  Eq.  (22)  leads  to 

A(m,n)  =  aMN8(m)8(n)  +  bM 


N- 1 

E  le 
/=  o 


— 27t nl  /N 


5{m). 


(30) 


The  ranges  of  indexes  in  Eqs.  (29)  and  (30)  are  m  =  0, 1, . . .  ,M  —  1  and  n  —  0, 1, ...  ,1V  1. 


Figure  3  is  the  matrix  structure  of  the  transform  in  Eq.  (29).  The  transform  is  of  a  very  special 
form  where  nonzero  entries  appear  only  in  the  first  column.  These  are  given  by 


A  (m ,  0)  =  aN 


A  (0, 0)  =  aN- 


MjM_  -  1) 

2 


M  ,  .M 
—  —  + 1—  cot 
2  2 


+  bMN 

(31) 

(m  =  1,2,. 

. .  ,M  —  1). 

(32) 
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Fig.  3  —  Matrix  form  of  the  two-dimensional  transform  A(m,n)  of  the 
function  X (kj )  =  ak  -f  h  showing  its  structural  aspects.  The  elements 
of  the  first  column  are  given  in  Eqs.  (31)  and  (32). 


Consequently,  even  when  M  —  N  we  obtain  a  very  asymmetric  transform  with  variable  entries  in  the 
first  column. 


The  matrix  structure  of  the  transform  in  Eq.  (30)  is  similar  to  that  of  Eq.  (29),  except  that  now 
we  have  the  nonzero  entries  appearing  in  the  first  row.  Figure  4  shows  the  structure.  The  elements 
appearing  in  the  first  row  are  given  by 


A  (0, 0)  =  aMN  +  bM  ^  ,  and 


A(0,n)  -  bM 


N  ,  .N 
—  + 1— cot 
2  2 


r  *>» 

7 xn 

) 

in  =  1,2,.  ..,2V  —  1) . 


(33) 

(34) 


It  is  clear  from  the  results  in  Eqs.  (31)  to  (34)  that  the  constant  terms  appearing  in  the  input 
functions  of  Eqs.  (21)  and  (22)  only  contribute  to  A  (0,0)  in  each  case.  They  only  offer  a  limited 
degree  of  freedom  so  far  as  changing  the  transform  and  need  not  be  included.  The  main  reason  for 
their  inclusion  here  is  to  demonstrate  what  the  general  linear  function  leads  to  under  the  discrete 
Fourier  transform  relation. 

CONCLUSIONS 

The  development  of  a  validation  procedure  for  the  FFT  algorithm  expresses  a  general  attitude  of 
the  author  that  algorithms  should  be  checked  by  application  to  a  case  for  which  the  results  are  known 
by  other  means.  This  is  not  always  done.  In  many  cases  the  reason  for  inadequate  tests  of  an  algo¬ 
rithm  may  be  traced  to  a  lack  of  known  analytical  results. 

We  have  presented  some  explicit  analytical  results  for  discrete  Fourier  transforms  in  the  one¬ 
dimensional  and  two-dimensional  situations.  The  results  are  easily  generalized  to  any  multidimen¬ 
sional  situation.  The  primary  motivation  for  the  work  was  to  establish  a  means  for  validating  multidi¬ 
mensional  FFT  algorithms  by  establishing  analytical  results  that  can  be  compared  to  the  output  of 
FFT  algorithms.  The  need  for  validating  FFT  algorithms  arises  in  many  contexts,  such  as  when  such 
an  algorithm  is  first  implemented  on  a  computer  facility  or  when  modifications  are  made  to  take 
advantage  of  new  techniques  or  new  devices  for  the  calculation  of  the  discrete  Fourier  transform. 


8 


NRL  REPORT  9042 


Fig.  4  —  Matrix  form  of  the  two-dimensional  transform  A(m,n)  of  the 
function  X(kJ)  =  a  +  bl  showing  the  structural  aspects  of  the 
transform.  Elements  appearing  in  the  first  row  are  given  by  Eqs.  (33) 
and  (34). 


The  general  idea  discussed  in  this  report  is  that  of  using  simple  discrete  functions  whose  discrete 
Fourier  transforms  can  be  evaluated  by  direct  analytical  means.  Several  functions  were  considered 
explicitly  for  the  one-dimensional  and  two-dimensional  cases,  and  their  advantages  and  disadvantages 
were  examined  in  terms  of  a  validation  procedure.  We  note  that  such  a  procedure  allows  the  algo¬ 
rithm  to  be  viewed  as  a  black  box  with  only  the  input  and  output  being  of  primary  interest. 


The  types  of  functions  we  have  dealt  with  are  of  the  form 

X(k)  =  ak  +b  ( k  =  0, 1,.. .  ,Af  -  l)and 


X(k,l)  —  akl  T* bk  +  cl  -\-d 


k  =0, 1, . . . ,  A/  —  1 

/  =  0, 1, . . .  —  1  J  ’ 


where  one  or  another  of  the  constants  may  be  zero.  These  functional  forms  appear  to  offer  suffi¬ 
ciently  general  tests  of  FFT  algorithms  in  one  and  two  dimensions. 


One  could  generalize  the  procedure  by  dealing  with  functions  of  the  form 


X(k)  »/(* )  (k 
X(k,  l)  =  f  (k)g(l) 


=  0,1,..., M-l)  and 

" k  =  0,1,...  ,M  -1 
/  =  0,1,...  ,N  -  1 

V  J 


and  so  on  for  higher  dimensions.  A  general  property  of  such  functions  is  that  the  multidimensional 
transform  of  such  separable  functions  will  be  found  to  reduce  to  a  product  of  one-dimensional 
transforms.  One  would  then  select  f(k),g(l)  so  that  their  one-dimensional  transforms  can  be  ex¬ 
plicitly  evaluated  in  analytical  terms.  The  only  difficulty  in  carrying  out  such  a  procedure  is  the 
analytical  evaluation  of  the  one-dimensional  transforms. 
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Appendix 

ANALYTICAL  EVALUATION  OF  A  DISCRETE  FOURIER  TRANSFORM 


We  here  consider  the  analytical  evaluation  of  a  discrete  Fourier  transform  of  a  specified  function 
X(k )  for  k  =0, —  1.  Despite  the  treatment  of  the  general  theoretical  aspects  of  such 
transforms,  the  evaluation  of  specific  results  is  rarely  treated  and  there  appears  to  be  no  extensive 
tabulation  of  transforms  of  simple  functions  that  can  be  worked  out  analytically.  The  particular  func¬ 
tion  we  are  interested  in  is  given  by 

X(k )  =  k  (k  =  0,1,...  ,M  -  1). 


It  is  desired  to  evaluate  the  transform  that  can  be  written  as 

M- 1 

A{m)  =  £  ke-lrikm,M  ( m  =  0, 1 ,...  ,M  -  1). 

k=  0 


(Al) 


For  the  sake  of  later  convenience,  we  note  the  summation  in  Eq.  (Al)  can  be  easily  evaluated  for  the 
special  case  m  =  0  since 


k  =0  Z 


(A2) 


Thus  we  now  have  to  evaluate  the  expression  in  Eq.  (Al)  for  m  =1,2, ...  ,M  —  1. 


To  carry  out  the  summation  in  Eq.  (Al),  it  is  useful  to  regard  the  discrete  functions  in  the  sum¬ 
mand  as  defined  for  all  values  of  k .  Thus  for  each  value  of  m  we  have  a  summation  problem  where 
it  is  desired  to  obtain  a  simple,  anlaytical  expression  for  the  sum.  Define  the  finite  difference  opera¬ 
tor  A  by 

A fi(k)  =  fi(k  +  1)  —  n(k) , 

and  note  that  it  has  the  property 

A  (jiv)  =  vAv  +  i >(k  +  1)A  n ,  (A3) 

where  it  is  understood  the  independent  variable  is  k  in  those  terms  where  it  is  not  indicated  explicitly. 
We  specify  that 

H(k)  =  k  ,  (A4) 


and 

A  i>  (k)  =  e~2*ikm/M  (k  =  0,1,..., M). 

This  allows  the  summation  problem  in  Eq.  (Al)  to  be  expressed  as 

m- t 

A(m)  =  £  n(k)  Av , 

k  =0 

which,  by  use  of  Eq.  (A3),  can  be  written  in  the  form 

M- 1  M  —  1 

A(m)  =  £  A  (nv)  -  £  v  (k  +  1), 

k  =0  k=  0 

or,  since  the  first  summation  can  be  summed  explicitly,  we  have 

m  — t 

A (m)  =  n{M)  v(M)  —  n(0)  v(0)  —  ^  +  1) . 

k  =0 


(A5) 


(A6) 
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The  reason  for  the  change  in  summation  is  that  the  summation  problem  in  Eq.  (A6)  will  be  able  to  be 
done  in  explicit  terms.  The  method  being  used  is  analogous  to  integration  by  parts  for  continuous 
functions. 


We  now  determine  v(k)  which  satisfies  Eq.  (A5).  We  can  define 

v(0)  =  0 

since  Eq.  (A5)  is  only  a  requirement  on  differences  of  v(k). 

Then,  we  find  by  evaluating  the  first  few  terms  of  the  expression  in  Eq.  (A5)  that 

„(1  )=e-2rikmlM  |fc=()> 


(A7) 


and 


_  -27 rikm/M  I  .  ^  —lirikm/M  \ 

v{Z)  -  e  \k=0  +  e  \k  =  l- 


On  the  basis  of  these  results  we  can  guess  that 

fc  —  i 

v(k)  =  J2  e~2nml/M(k  =  1,2,...  ,M). 
/  =  0 


(A8) 


This  can  be  verified  by  forming 

which  gives 

and 


A  v(k)  =  vQc  +  1)  -  v(k), 
A  v  (0)  =  ^(1)  =  e~2iriml/M 


i=o > 


(A9) 


Jc  k  —  1 

y(k  +  1)-  v(k)  =  fj  e~2riml/M  _  £e-2 *iml/M 
l ~0  1=0 


or 


A  v  (k)  =  v  (k  +  1)  —  v  (k)  =  e 


-lirikm  /M 


(k  =  1,2 


(A10) 


It  is  clear  that  the  functions  defined  by  Eqs.  (A7)  and  (A8)  solve  the  difference  equation  exhibited  in 
Eq.  (A5). 

We  now  find  an  explicit  expression  for  the  summation  in  Eq.  (A8).  Write  Eq.  (A8)  as 
k- 1 

v(k )  =  (cos(27rm//M)  —  /  sin(27rm//M))  (k  =  1,2  (All) 

/  =o 

The  summations  appearing  in  Eq.  (All)  are  known  (Hamming  1962,  p.  44)  and  m  ^  0,  can  be  writ¬ 
ten  as 


k- 1 

J2  cos  (27rm/  /M)  = 
/  =o 


.  27rm 

r  *> 

,  i  ,  1 

+  sin 

r- 

Tim 

sin 

k  —  1  H - 

M 

2 

j 

M 

2  sin 


7r  m 


M 


(A  12) 


and 


k  —  1 

^2  sin(27rm//M) 
/=o 


cos 


2irm 

M 


k  —  1  H - 

2 


+  cos 


urn 


2  sin 


7 xm 


(A  13) 


12 


NRL  REPORT  9042 


where  k  =  1,2 ,M .  This  completes  the' explicit  determination  of  the  function  v(k). 


We  now  return  to  Eq.  (A6)  to  find  what  simplifications  can  be  made.  We  note  that  /t(0)  and 
p(0)  both  vanish.  Also,  v{M)  can  be  evaluated  from  Eqs.  (All)  to  (A13).  We  have 


.  27 xm 
Sm  M 


u(M)  = 


r  "> 

,  i 

c 

Tim 

2  Tim 

,,  ,  1 

M  -  1  +  — 

+  sin 

+  i  cos 

M  -1  +  — 

2 

M 

M 

2 

-J 

- 1  cos • 


7 vn 
M 


2  sin 


7rm 


and  on  simplification  of  the  trigonometric  terms,  the  result  is  obtained  that 

v(M)  =  0. 

Consequently,  Eq.  (A6)  takes  the  form 

M- 1 

A(m)  =  -  £  v(k  + 1) 


k  =0 


M- 1 

■  £ 

k=  0 


.  2ir  m 

sin  ~TT 
M 


k  + 


1 


+  sin  ■ 


Tim 

2  Tim 

i  ,  i 

Tim 

M 

C0S  HM 

M 

k  2 

—  cos 

M 

_L  i 

L  J 

L  J 

2  sin 


7 im 

M 


2  sin 


7rm 

M 


•  (A  14) 


We  now  consider  each  term  in  Eq.  (A14)  separately.  We  obtain 


1  .  2-Ktn 
£  sin- 

k  =0 


M 


k  + 


cos 


2  Tim 

M 


M-  1+-^ 
2 


+ 


7Tff7 


cos 


7r  m  7r  m 

M  M 


2  sin 


7r  m 


which  is  found  to  vanish  on  simplification.  Likewise,  we  get 


2-irm 

£  cos- 
k  =0 


M 


k+2 


=  0. 


Eq.  (A14)  then  simplifies  to 


or 


M- 1 

AQn)  =  -  £ 

k  =0 


..  x  M  ,  .  M 
A(m)  =  -  —  +  i  —  cot 


“ 

1 

.  1 

Tim 

— 

—  1  —  cos 

— — — 

2 

2 

M 

J 

nm 


M 


(m  =  1,2,. . .  ,M  —  1). 


This  is  to  be  supplemented  by  the  results  indicated  in  Eq.  (A2)  that 

A(0)  =  mLzA' 


(A15) 


(A16) 


The  expressions  given  in  Eqs.  (A  15)  and  (A  16)  provide  the  explicit,  analytical  expression  for 
the  discrete  Fourier  transform  of  the  function  X(k)  -  k,(k  =  0, 1,. . .  ,M  -  1),  which  is  the  result 
desired.  It  should  be  clear  that  even  relatively  simple  discrete  functions  lead  to  rather  difficult  sum¬ 
mation  problems  when  it  is  desired  to  obtain  an  explicit  representation. 
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Figure  A1  shows  a  plot  of  the  imaginary  part  of  the  discrete  Fourier  transform.  When 
normalized  by  M  the  transform  values  fall  on  a  continuous  curve  that  is  independent  of  M.  The  real 
part  of  the  transform  is  of  a  simple  form  and  is  not  exhibited  in  graphical  form. 
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M 


Fig.  A1  —  The  imaginary  part  of  the  one-dimensional 
discrete  Fourier  transform  A(m)  of  the  discrete  function 
X(k)  =  k(k  =  0, 1,*  *  *  ,Af  —  1).  It  is  plotted  in  the 


form  Im 


iA{m) 


because  when  normalized  by  M 


the  discrete  values  fall  on  a  single  continuous  curve  that 
is  independent  of  M.  The  solid  circles  are  discrete 
values  of  the  imaginary  part  of  the  normalized  transform 
for  the  case  when  M  —  2n ,  where  n  >  4.  For  other 
choices  of  M ,  the  values  would  still  be  along  the 
continuous  curve  but  not  always  at  the  circled  points. 
The  imaginary  part  of  the  transform  always  vanishes  at 

m  =  0.  The  continuous  curve  is  given  by  cot 

M 

j 

with  m  regarded  as  continuous.  Values  of  the  imaginary 
part  of  the  normalized  transform  outside  the  domain  can 
be  obtained  by  periodic  extension  of  the  values  shown. 
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